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Abstract. We consider the numerical evaluation of the Evans function, a 
Wronskian-like determinant that arises in the study of the stability of travel- 
ling waves. Constructing the Evans function involves matching the solutions 
of a linear ordinary differential equation depending on the spectral parameter. 
The problem becomes stiff as the spectral parameter grows. Consequently, the 
Gauss— Legendre method has previously been used for such problems; however 
more recently, methods based on the Magnus expansion have been proposed. 
Here we extensively examine the stiff regime for a general scalar Schrodingcr 
operator. We show that although the fourth-order Magnus method suffers from 
order reduction, a fortunate cancellation when computing the Evans matching 
function means that fourth-order convergence in the end result is preserved. 
The Gauss-Legendre method does not suffer from order reduction, but it does 
not experience the cancellation either, and thus it has the same order of conver- 
gence in the end result. Finally we discuss the relative merits of both methods 
as spectral tools. 



1. Introduction 

Many partial differential equations admit travelling wave solutions; these are 
solutions that move at a constant speed without changing their shape. Such trav- 
elling waves occur in many fields, including biology, chemistry, fluid dynamics, and 
optics. It is often important to know whether a given travelling wave is stable: 
does it persist under small perturbations? A major step towards determining the 
stability of a travelling wave is to locate the spectrum of the linearization of the 
differential operator about the travelling wave. Evans [TU] considers a shooting 
and matching method for this task. Evans introduced a function to measure the 
mismatch for a specific class of reaction-diffusion equations. This function was 
called the Evans function by Alexander, Gardner and Jones [2], who generalized 
its definition considerably. Since then, the Evans function has been used frequently 
for stability analysis; see, for example, [TJ [3J HI Q~5] for numerical computations 
employing the Evans function and [Til H21 [30l [34l [35] for an analytic treatment. 
The review paper by Sandstede 33J gives an excellent overview of the field. 

The Evans function is a function of one argument, the spectral parameter A, 
and zeros of the Evans function correspond to eigenvalues of the corresponding 
operator. Hence, one can get information about the spectrum by finding zeros of 
the Evans function, either analytically or numerically. Our focus here is on the 
numerical approach. 
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The main part in the numerical evaluation of the Evans function is the solution 
of a linear ordinary differential equation depending on the spectral parameter A. 
This is often done with an off-the-shelf integrator using an explicit Runge-Kutta 
method. Afendikov and Bridges [I] noticed that when A grows, the problem may 
become stiff, and therefore they use a Gauss-Legendre method. Recently Aparicio, 
Malham and Oliver 3 proposed a new procedure based on the Magnus expansion, 
building on the work of Moan [25] and Greenberg and Marietta [14] who used the 
Magnus expansion to solve Sturm-Liouville problems. Aparicio et al. noticed that 
the Magnus method suffers from order reduction in the stiff regime. This means 
that the fourth-order integrators whose global error should scale like h 4 when the 
step size h is small, instead converge more slowly. They analyzed this phenomenon 
in a modified Airy equation using the WKB-method. However, they restricted 
themselves to those values of A which correspond to the essential spectrum of the 
linearized differential operator. 

The current paper continues the analysis of the Magnus method in the context 
of Evans function evaluations. We concentrate on scalar Schrodinger operators to 
simplify the analysis. Other methods based on a transformation to Priifer vari- 
ables [32J probably perform better in the scalar setting, but these methods cannot 
be used unchanged in the non-self-adjoint case where our interest lies. 

We present another approach to the analysis of the Magnus method based on a 
power series expansion, which is valid for values of A outside the essential spectrum. 
We will show that the Magnus method also suffers from order reduction in this 
regime. Specifically, the relative local error is of order A _1 / 2 /i 2 as h ~ * with 
|A| 1//2 ft. 1. However, there are two subsequent important observations. Firstly, 
when going from the local to the global error, one does not lose a factor of h (as 
usually), but the global error is also of order A _1 / 2 /i 2 . Secondly, the order reduction 
disappears completely when we evaluate the matching condition: the relative error 
in the Evans function is of order A _1 / 2 /i 4 , thus quartic in the step size, just as one 
would expect from a fourth-order method. Since useful asymptotic estimates invoke 
an order A -1 error, at best, our numerical schemes even with order reduction prove 
a useful spectral tool in the regime |A| <C h~ s . 

The phenomenon of order reduction was discovered for implicit Runge-Kutta 
methods by Prothero and Robinson 31J. Nowadays, it is understood within the 
framework of B-convergence (see for instance jTHl §IV.f5]). The stability of Mag- 
nus methods has been analyzed for highly-oscillatory equation by fserles [19], for 
Schrodinger equations by Hochbruck and Lubich [17] . and for parabolic equations 
by Gonzalez, Ostermann and Thalhammer [13] . Unfortunately, these results cannot 
yet be fitted into a general theory [5D] . The present paper can also be viewed as a 
contribution to this research. 

We also present an analysis of the fourth-order Gauss-Legendre method. We 
show that the relative error committed by this method does not contain a term 
of order A _1 / 2 /i 2 , but only smaller terms. Furthermore, the error decreases even 
further when evaluating the Evans function. As explained in more detail later, this 
is due to an effect similar to the one which makes the trapezoidal rule very efficient 
for the quadrature of periodic functions. 

The contents of this paper are as follows. In the next section, we define the 
Evans function and we give an asymptotic expression for the Evans function in 
the scalar case when the spectral parameter A is large in modulus and outside the 
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essential spectrum. We then define the Magnus method in Section [3l We show 
that the Magnus method suffers from order reduction and we compute the error 
when evaluating the Evans function. We repeat the computation for the Gauss- 
Legendre method in the next section. The analysis is corraborated by numerical 
experiments in Section [5] In the final section, we compare our results with those of 
Aparicio, Malham and Oliver [3] and we discuss the stability of the Magnus method 
in general. More details of the intricate calculations presented in Sections [2HH can 
be found in the technical report |29j . 

2. The Evans function 

We are interested in homogeneous reaction-diffusion equations on an unbounded 
one-dimensional domain. Such equations have the form 

(1) u t = Ku xx + f(u), 

where K is an n-by-n diagonal matrix with positive entries (the diffusion coeffi- 
cients) and the unknown u is a function of t and x. The function / : R™ — > R" 
describes the reaction term; we assume that / is sufficiently smooth. 

A travelling wave solution has the form u(x,t) = w(£) with £ = x — ct where 
c is the wave speed — see for example Kolmogorov, Petrovsky and Piskunov [23j . 
We assume that a travelling wave solution for the equation is known, at least 
numerically. Furthermore, we assume that u is constant at infinity, meaning that 
the limits u± = \im^± 00 u^) exist (in fact, we will need later that additionally, 
the derivatives flW vanish at infinity for p = 1,2,...). Such a wave is called a pulse 
(if u + = U-) or a front (if u + ^ u_). 

To study the stability of the travelling wave, we linearize |T]) about the wave and 
write the result in the (£, t) coordinate system which moves with the same speed 
as the travelling wave. This yields 

(2) u t = Ku££ + cu^ + Df(u) u. 

Define the operator C by C{U) = KU" + cU' + Df{u) U, where £/ : R — > C™. Its 
spectrum determines whether @ has solutions of the form u(£,i) = e xt U(£). The 
stability of the travelling wave u can be deduced from the location of the spectrum 



The spectrum of C can be divided in two parts: the point spectrum cr pt (£), 
consisting of those A € er(£) for which C — XI is Fredholm of index zeroQ and the 
essential spectrum <7 CS s(£), which contains the rest of the spectrum. We assume 
that <7 ess (£) is contained in the left half- plane {z £ C : Re z < 0}. This means that 
the spectral stability is determined by the position of the eigenvalues A € a pt (£). 

We now introduce the Evans function, which is a tool for locating these eigenval- 
ues. We rewrite the eigenvalue equation C(U) — XU as the first-order differential 
equation 



of C. 



(3a) 




An operator is Fredholm with index zero if its range is closed and the dimension of the null 
space equals the codimension of the range. 
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where y : R — ► C 2 ™ and the matrix A is given by 



4(6 A) - [ /f -i( A/ _^ /( y 



(0)) -cK- 1 



Since the spectral problem (|3a|) is a linear equation, its solutions form a linear 
space of dimension 2n. Define E~(X) to be the subspace of solutions y satisfying 
the boundary condition ?/(£) — > as £ — > — oo. Similarly -E + (A) denotes the 
subspace with y(£) — > as £ — + +oo. Any eigenfunction must satisfy both boundary 
conditions and hence lie in the intersection of E~(X) and E + (X). 
For all A ^ cr ess (£), we have 



Choose a basis yi( ■ ; A), . . . ,yk( ■ ; A) of £7~(A), where = dimi? _ (A), and a basis 
y/s+i( ■ ; A), . . . , y2n( ■ ; A) of E + (X). We can assemble these basis vectors, evaluated 
at an arbitrary point, say £ = 0, in the 2n-by-2n matrix 



The Evans function, denoted D(X), is defined to be the determinant of this matrix. 
If the determinant vanishes, then the yi are linearly dependent, which implies that 
the spaces E~{X) and E + {X) have a nontrivial intersection, and this intersection 
contains the eigenfunctions of (f3aj) . Therefore, 13(A) = if and only if A G a v %(£). 

Let C denote the connected component of C \ er css (£) containing the right half- 
plane. We can choose the basis vectors yi to be analytic functions of A in the 
region C. The Evans function will then also be analytic in C and the order of its 
zeros corresponds with the multiplicity of the eigenvalues of C 

More details on the Evans function and the stability of travelling waves can be 
found in the landmark paper by Alexander, Gardner and Jones [2] and the review 
article by Sandstede [55] , 

2.1. The Evans function near infinity. We are interested in the behaviour 
of D(X) and numerical approximations to D(X) as |A| — * oo, because experiments 
show an unexpected deterioration of the approximations in this limit [3] . For sim- 
plicity, we will restrict ourselves to scalar reaction-diffusion equations, i.e., we as- 
sume that n = 1. However, it is expected that the methods of analysis presented in 
this paper also apply to the nonscalar case, though the computations will obviously 
be more involved. 

We may assume without loss of generality that the diffusion coefficient is 1, so 
that the partial differential equation reads 



dim£T(A) + dimE + (X) = 2n. 



yi(0;A) ... y k (0-X) y fc+ i(0;A) ... y 2 „(0;A) 




where 



(4b) 







1 



A 



— c 
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±00 are given by 

1 ' 

A-/'(fi±) -c 



The limits of A as £ 

(5) A ± (X) = 
Furthermore, the eigenvalues of A_(A) are 

(6) M L 1] , mL 2] = i (-c± Vc 2 + 4(A-/'(,i_))) . 

To avoid any confusion between the eigenvalues of the differential operator C, which 
form the point spectrum that we want to compute, and the eigenvalues of the 
matrices A±(X), we call the latter spatial eigenvalues. 

One of the spatial eigenvalues is purely imaginary if A lies on the parabolic curve 
given by 

7 _ = { -s 2 + f'(u_) + ics : s e R }. 
The curve 7_ is part of the essential spectrum. If A lies to the right of t_, then the 
spatial eigenvalues /J?} and have positive and negative real parts, respectively. 

The limit £ — > +00 is treated in the same manner and leads to the curve j + . 
The region C on which the Evans function is defined is the part of the complex 
plane to the right of 7_ U 7+ . 

We assume henceforth that A e C. We compute the asymptotic behaviour of the 
Evans function as |A| — > 00 using a different approach to that outlined in Alexander, 
Gardner and Jones [2, §5B], extending the approximation to further higher order 
corrections. We start with the solution y of ^ satisfying $/(£) — * as £ — > —00. 
The matrix A in (|4b[) goes to A_ as defined in ([5]) in this limit, and the eigenvalues 
of A- are given in 
This suggests writing y as 



with corresponding eigenvectors (1, /jP}) t and (1,/i^ 2 ') 



(7a) 2,(0 = cxp( M W £)(n(£) 
where 
(7b) 



1 



1 

[2] 
I*- 



= exp(^ ] OBj/(0 



It 




1 1 " 


and £> 


[1] [2] 


u 







The vector y satisfies the linear differential equation 
(8) = i(£;A)y with A={B~ 1 AB- 

The matrix j4(£; A) in this equation is given by 

-£M0 -^-(0 



(9a) 

where 
(9b) 



MO = /'(fi(0) - /'(«-) and « = Vc 2 +4(A -/'(«_)). 



Note that the parameters c and A are replaced by only one parameter, k. Now, 
suppose that u and v can be expanded in inverse powers of k: 
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If we substitute these expansions in ([8]) and equate the coefficients of the powers 
of K, we find: 



= 0, 
i' = 0, 



= 



-Vq, 

-n, 



u[ = -tp-(g) (uo + «o), v'i 
u' 2 = -ip-(£) (wi +vi), v' 2 

Assuming that and v(£) are bounded as £ 
tions (up to a multiplicative constant) is 



(10) 
where 



u(£;/s) = 1 - + 

i>feK) = K-V(0+O(K- 3 ). 



<P-(0 ("0 + «o) - V 2 , 

— oo, the solution of these equa- 

(<M0) a + o(«- 3 ), 



M0 



(a;) dec. 



We can do something similar to find the solution y of @ satisfying ?/(£) — > as 
£ — > +oo. Instead of (JT)), we write y as 



2/(£) = exp(M+0 •£+!/+(£) where R, 
Expanding y + in negative powers of where 



[2] 



(11) k+ = v/c2+4(A -/'(«+)), 

similar to (I9bl), we find that 



(12) y+ = 
where 



with 



'fi+K;/s) = K7V+(0 + O(Kl 3 ) 



+ h 

v+(£; k) = 1 + kT^MO + K 2 (<M£)f + 0(k- 3 ), 



¥'+(0 = /'(fi(0)-/'(fi+) and MO 



¥>+(a5) dx. 



The Evans function is obtained by evaluating both the solution satisfying y(£) — > 
as £ — * — oo and the one satisfying y(£) — ► as £ — > +oo at £ = 0, collecting the 
resulting vectors in a matrix and computing the determinant of this matrix. This 
yields 



D(X) = (By(0)) A {B+y+(0)) 



(13) 



ti(oy 

5(0) 



\(k — c) — \{n + c) 
= i(«-K + )(tj(0)C+(0)-fi(0)«+(0)) 

+ |(« + «+)(i;(0)fi+(0)-u(0)t; + (0)) 



1 

) -i(K+ + C ) 



't»+(oy 

«+(0) 



Substituting (fT0|) and (fT2|) and using the fact that k—k + — 0(|A| ^ 2 ) as |A| — > oo, 
we find that 



(14) 



D(A) = K+)i2(0)u + (0) + 0(kT 2 ) 



= -2A 1/2 + $ - iA" 1 / 2 ^ 2 - 2/'(fi_) - 2f'(u+) + c 2 ) + ^(A- 1 ), 
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where 



3> = $_(0) +$+(0) 



(f- (x) dx ■ 



ip+(x) dx. 



The approach of Sandstede [33] yields D{\) = -2A 1 / 2 +C(l). This agrees with flU]) , 
but the approach presented here gives two more terms. Furthermore, we can easily 
find additional terms by extending the expansions (jlpp . 



3. Magnus methods 

If we want to evaluate the Evans function numerically, we have to solve the 
differential equation ([3J . Moan [25] studied methods based on the Magnus series for 
the solution of Sturm-Liouville problems of the form — (py 1 ) 1 + qy — \wy on a finite 
interval. Moan noticed that some of the quantities involved in the computation are 
independent on the spectral parameter A and hence need to be computed only once 
when solving the differential equation for several values of A. Moan also proposed a 
modification of the method based on summing some of the terms analytically which 
improves the accuracy when A is large. 

Jodar and Marietta [5T] noticed that the Magnus method in combination with the 
compound matrix method performs well on some scalar Sturm-Liouville problems 
of high order; see also Greenberg and Marietta [14] . 

This approach was generalized by Aparicio, Malham and Oliver [3J, who pro- 
posed to use a Magnus method for solving the boundary value problem (J3J. They 
mentioned the robustness across all regimes as an advantage of Magnus integrators. 
Furthermore, they pointed out that the computational cost of Magnus methods, as 
well as other methods, can be decreased by means of a precomputation technique, 
In this section, we further analyze the behaviour of the Magnus method in the 
regime where A is large in modulus. 

Magnus [24] showed that the solution of the differential equation y' = A(£)y can 
be written as = exp(f2(£)) y(0), where the matrix i s given by the infinite 



A(x) dx 



(15) 



? r 





xi 



A(x 2 )dx 2 ,A(x 1 ) 



dxi 



A(x 2 ) dx 2 , 



A(x 2 )dx 2 , A{x\) 



A(x 3 )dx 3 , A(x 2 ) 



dx 2 , A(xx) 



dx\ 



dx\ 



where [ • , • ] denotes the matrix commutator defined by [X, Y] = XY — YX. Moan 
and Niesen [26] proved that the series converges if J Q ||^4(x)|| dx < ir. 

The Magnus series can be used to solve linear differential equations numerically, 
if we truncate the infinite series and approximate the integrals numerically. For 
instance, if we retain only the first term in the series and approximate A(x) by the 
value at the midpoint, we get f2(£) ~ hA(^). The resulting one-step method is 
defined by 



(16) 



y k+ i = exp(hA(£ k + \h)) y k , 



where h denotes the step size and y k approximates the solution at = £o + kh. 
This method is called the Lie midpoint or exponential midpoint method. It is a 
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second-order method: the difference between the numerical and the exact solution 
at a fixed point £ is 0(h 2 ). 

We can get a fourth-order method by truncating the Magnus series (TT5|) after 
the second term. We replace the matrix A(£) by the linear function Ao + ^Ai which 
agrees with A(£) at the two Gauss-Legendre points 

(17) $ ] = & + (| - ±VS)h and £[ 2] = & + (| + iV3)A. 

This yields the scheme 

(18a) y fc+ i = exp(n k )y k , 

where 

(18b) = Jfc(A«fl) + A(ef")) - ^ 9 [A(fW),A(eP")]. 

The reader is referred to the review paper by Iserles, Munthe-Kaas, N0rsett and 
Zanna |20j for more information on Magnus and related methods. 

If we define yk by y k = exp(//[^£) By k , as suggested by (J7J, then the fourth-order 
method ([TBI) transforms to 



(19a) y k +i = exp(n fe ) y k , 

where 

(19b) ^ fc = i h (A^) + A($)) gh'[A{tP),A(t]P)], 

with A as given in @. So applying the Magnus method to the transformed equa- 
tion {8| and transforming the result back to the original coordinate system is the 
same as applying it to the original equation. This can be explained by the equivari- 
ance of the Magnus method under linear transformations and exponential rescal- 
ings [9]. 

Substitution of © in (|T9b"| yields 



-K 1 a k (3 k - k 1 a k 
f3 k + K~ 1 a k -K + K,- 1 a k 



(20a) n k = h 

with 

(20b) a k = fo4$ ] )+<p-($)) and & = -^h{^[ 1] ) - 

Note that a k and (3 k approximate if- and j^h 2 (p'_ respectively. Furthermore, the 
exponential midpoint rule (fT6|) is also given by (|20a| . but with a k — (f-(£k + 2^) 
and (3 k = instead of (|20b|) . 

3.1. Estimates for the local error. The local error of a one-step method is the 
difference between the numerical solution and the exact solution after one step. For 
the Magnus method, the local error is 

L k = cxp(n k ) y(£ k ) - y(£fc+i), 

or, in transformed coordinates, 

(21) L k = exp(^fc) - y(£ k+ i). 

The exponential of the matrix fl k is most easily calculated by diagonalization: if 
Qfe = VfcAfcVjT 1 with diagonal, then exp(fifc) = V k exp(A^) V k ~ x and exp(Afc) is 
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formed by simply exponentiating the entries on the diagonal. In fact, the diagonal 
entries of Ak are 



(22) 



as \K\ 



=h[ K - 1 ((3 2 k -a k ) 



A" = -h 



(Pi- 



ak) - k 



a k 



5 (A 2 



0(k 



-5- 



ak 



0(k 



The definition of 



9b|) implies that Re n > unless A is real 



and A < /'(tt_) — jc 2 . Under this condition, — RcAjJ ^ 1 if h\n\ ^ 1 and hence 

exp(A). ) is exponentially small. 

We now assume that we are in the regime with |A| 3> h~ 2 , h — > 0, and A bounded 
away from the negative real axis in the sense that | arg A| < n — e where e > 0. In 

[2] 

this regime, exp(A^ ) is exponentially small. Taking this into account, a lengthy 
but straightforward calculation shows that 



(23) exp(ft fc ) 



1 



K 



Pk 
K 



I3 k a k + hf3kXk 



ak + hfi k Xk 
hkxk , 



0(«- 4 ) 



where Xk = a k — P 2 - Substituting this result and the approximation (|10[) for the 
exact solution in the definition (|2"Tj) . and using the definitions of Xk, Pk, and <E>_, 
we find that the local error of the Magnus method is given by 

'n-^k + Oi^h 4 ) 

K~ X p k + 0(K- 2 h) 



(24) 



U 



where 
(25) 

and 
(26) 



Ik = 



ip-(x) dx - h(a k - Pi) 



h5 (^o^'Kfc + \h) + M^'ttk + \h)) 2 ) + o(h 7 ) 



Pk = ±h 2 v '{tk + \h) + 0(h 4 ). 

In deriving the above expression, we also replaced by ip(£) = f'(u(^)). This 

is allowed since they differ by a constant term, and only the derivative appears 
in (EU). 

Equation gives the local error of the fourth-order Magnus method (|T5)l in 
transformed variables. Since the method has order four, the local error is 0(h 5 ) as 
h — > when solving a fixed equation. However, in our case, the constraint | A| 3> hr 2 
implies that A, and hence the relative influence of the coefficients in the equation, 
must change as h approaches zero. It turns out that the local error is 0(h 2 ) in this 
setting. In other words, the method behaves like a first-order method (globally). 
This phenomenon is called order reduction. 

The cause of this order reduction is the stiffness of the differential equation J4} . 
Indeed, if we define the stiffness ratio as the quotient between the largest and 
smallest eigenvalue (as in Iserles [IS]), then the stiffness quotient is 



•0(1), 



aL 2] 




K 2 


aL 11 




Pi - a k 
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where the leading term shown grows like |A|. Hence, the problem is stiff if A is large 
in modulus, causing troubles for the numerical method. 



3.2. Estimates for the global error. The global error is the error of the numer- 
ical method after several steps, say k. Hence, the global error is E k — y k — 2/(£fc)j 
with y k defined by the numerical method starting from yo = 2/(£o)- For the Magnus 
method (Tig)) , the global error satisfies the recursion relation E k +i = cxp(Slfc) E k +L k 
with Eq = 0, or, in transformed coordinates, 



(27) 



E k+1 = en-p(Q k ) E k + L k , E o = 0. 



A routine induction argument using 
the global error is given by 



and 



shows that the leading term of 



(28) 



E k = 



k _1 A-i +<D{n- 2 h) 



This shows an advantage of stiffness: the exact flow quickly reduces the error in the 
stiff component. The Magnus method inherits this property here and annihilates 
at every step the error in the stiff component up to leading order (if \n\h 1). On 
the other hand, the first (nonstiff) component of the error is propagated without 
change. Since the local error in the stiff component is much bigger than the error 
in the nonstiff component (order h 2 versus order h 5 ), we arrive at the surprising 
conclusion that the local and global error are equal at leading order. 

Substituting the definitions of a k , (3 k , and ^y k in (|20b[) and (|2B[) and approximat- 
ing the sum by an integral, we find that 
(29) 



E k — 



We see that the global error is 0(h 2 ). Usually, a factor h is lost in the transition 
from the local to the global error, but here both the local and the global error 
are 0(h 2 ), because the local error is mainly in the stiff component. Hence, the 
fourth-order Magnus method given in (fT8|) behaves like a second-order method if 
one considers the global error. The numerical experiments in Section [5] support this 
analysis. However, the fact that the global error is 0(h 2 ) does not tell the whole 
story, as the relative error (the error divided by the magnitude of the solution) 
may give a better picture. Indeed, it follows from ([7]) and (fTOf that the solution y 
grows as \/A, so the relative error is 0(h 2 /V\). In other words, the relative error 
decreases as |A| — > oo. 

We can compute the error associated with the solution satisfying the boundary 
condition that y(£) — ► as £ — > +oo similarly. Instead of (|28|) . we now have 



(30) 



E i=0 r 



■0{n+ 2 h A ) 



where at , /3t and 7^ are as given in (|20b[) . with cp+ and <i> + replacing and $_ 
respectively, and k + is as given in (fTTj) . 
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3.3. The error in the Evans function. The Evans function given in (fT3|) is 

D(X) = (By(0))A(B + y+(0)). 
The numerical error when evaluating the Evans function is therefore 
(31) E D = (By(0)) A (B+E+) + (BE k ) A (B+y+(0)) + (BE k ) A {B + E+). 



We can expand this in the same manner as in (|13[) . The first term on the right-hand 
side becomes 

±(k - k + )([E] 2 v + (0) - [£]iti+(0)) + i( K + k+)([E] 2 u + (0) - [£]i«+(0)). 

The dominating term in this expression is — ^(k + k + )[E]iv + (Q), which is O(X h 4 ); 
all other terms are 0(X^ 1 h 2 ) or smaller (recall that we assumed that |A| ^ h~ 2 ). 
Therefore, the dominating contribution to the error in the Evans function comes 
from the first (nonstiff) component of the global error, even though the second 
(stiff) component is larger. This is because the stiff and nonstiff directions are 
exchanged when you integrate in the other direction. Hence, when taking the 
wedge product of the global error of the solution on [— oo,0] with the solution 
itself on [0, +oo], the stiff component of the global error is paired with the stiff 
component of the solution; similarly, the nonstiff component of the global error is 
paired with the nonstiff component of the solution. Since the solution is mainly 
along the nonstiff direction, the nonstiff component of the global error (which has 
order h ) is brought to the fore, and the stiff component of the global error (which 
has order h 2 ) is reduced. 

Substituting the exact solution from (|10[) and (jTSJ) and the global error from (l2"8)) 
and (j3"0|) in (|3~lj) . we find that the error in the Evans function is 

E D ^ 1J+ k ^i+o{\-" 2 h% 

3=0 j=0 

Assuming that the differential equation is solved on the intervals [— L, 0] and [0, L], 
with L — Nh, this evaluates to 



(32) 



JV-l fL 



j=-N 



fc-1 k-1 



The term within brackets is difference between the approximation of J_ L ip(x) dx 
by two-point Gauss-Legendre quadrature and the integral itself. In our setting, all 
derivatives of the travelling wave u, and therefore also of the function ip = f o u, 
vanish at infinity (see Now, assuming that L is so large that the derivatives 
of ip at L are negligible, the error in Gauss-Legendre quadrature vanishes at all 
orders in h, for essentially the same reason that the trapezoidal rule is so effective 
for periodic integrands; this is easily proved with the Euler-MacLaurin formula 
(see, for instance, Davis and Rabinowitz [71 §3.4]). Hence, only the sums involving 
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the (3j and (5^ survive. These can be approximated easily using (|26|) . and we find 
that 

(33) E D = -—J (<p'(x)) 2 dx + 0(h 6 ) with <p(t) = f'(u(0)- 

So, in the end, the error in the Evans function is of order ft 4 , which is just what 
one would expect from a fourth-order method. 

3.4. The exponential midpoint rule. We saw above that the fourth-order Mag- 
nus method (fT8|) suffers severe order reduction when solving (j4|) with |A| ^S> 1. This 
is not the case for all methods. There are even methods based on the Magnus 
series which do not suffer global order reduction, like the exponential midpoint 
rule (flU)) which has order two (this method is also known as the second-order Mag- 
nus method). When applied to this method is of the form p9aj) . (|20a[) with otk 
and Pk given by a k = <P-{£,k + 5^) an d At = respectively. If we substitute this 
in we see that the k^ 1 term in the second (stiff) component drops out, and 
that the local error is [C(k _1 /i 3 ) 0(k~ 2 /i)] t . So, the exponential midpoint rule 
does suffer some local order reduction, but not as severe as the fourth-order Magnus 
method, for which the second component of the local error is of order k _1 /i 2 . 

The global error can be computed as in £)3.2I Again, only the nonstiff component 
propagates, so the global error is [C>(k _1 /i 2 ) 0(K~ 2 h)] T . As \n\ ^> ft -1 , the first 
component dominates and the exponential midpoint rule effectively does not suffer 
from order reduction if one looks at the global error. 

Continuing to find the error in the Evans function, as we did in H3.31 we find 

N-l fL 

(34) E D = h v{jh+\h)- I <p{x) dx + <D{\- 1/2 h 2 ). 

j=-N J - L 

Comparing with (|3^|) for the fourth-order Magnus method, we see that the sums 
involving the (3j and (3^ have dropped out (because 0j = 0), and that the two-point 
Gauss-Legendre quadrature is replaced by the trapezoidal rule. Again, the error 
in the trapezoidal rule vanishes at all orders if L is sufficiently large. Hence, the 
error in the Evans function is C(A _1 / 2 /i 2 ). In contrast, the fourth-order Magnus 
method has Ed = 0(h 4 ), see (|3"3")l . Thus, we can expect the second-order method 
to be more accurate than the fourth-order method. The experiments in Section [5] 
confirm this. 

4. The Gauss-Legendre method 

Most numerical computations of the Evans function reported in the literature 
use a Runge-Kutta method, in particular the classical explicit fourth-order method 
and the two-stage Gauss-Legendre method. As the differential equation that we 
want to solve, is stiff, we consider the two-stage Gauss-Legendre method. The 
method is given by 

Sl = A(Z l k 1] )(y k + ±h Sl + (t-^)h S2 ), 

( 35 ) s 2 = A(£ ] ) (y k + ({ + %h Sl + ±hs 2 ), 

Vk+i = vk + lK s i + s 2), 

where £^ and ^ are the Gauss-Legendre points, given in (fT7|) . 
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We will now analyse the error committed by the Gauss-Legendre method when 
computing the Evans function. After the usual coordinate transformation, cf. (|7|). 
and substitution of the matrix given in @, we can solve the system After a 

lengthy but relatively straightforward calculation, we find that y k +i = ^kVk with 



1 — K h(Xk 



1 ^-2l2„,2 



1 



l2 K - 2 h- 1 p k 
YlK- x h- x + 72n- 2 h- 2 



0(k 



y-k i 2'" 

where a k and [3 k are as defined in (|20bj) . The local error can now be found by 
substituting this matrix and the exact solution in L k = ^ky(^k) — f/(£fc+l)i °f- (EH- 
This yields 



(36) 

where 
(37a) LI 

(37b) L\ 



L k = 



- 2 L b k +0(n- 3 h 



3 h 5\ 



y -L c k + 0(n- 3 h 2 ^ 



<f-{x) dx - ha k = 0(h 5 ) 



6, 



Uha k f 



tk+h ^ 2 

<P-(x) dx 



(37c) LI 



if-(x) dx 

-L a k (^ k ) + ha k + IL%) =0(h 5 ) 

12k- 1 k - p_ (& + h) + <p-{£ k ) = 0(h 3 ). 

For reasons which will soon become clear, we must retain the k~ 2 term in the 
above expression, in contrast to (|24|) for the Magnus method. We see that the two 
terms in the first (nonstiff) component are of order K~ x h 5 and K~ 2 h 5 , while the 
second (stiff) component is of order n~ 2 h 3 . This is a similar situation as with the 
exponential midpoint rule, except that (|35|) is a fourth-order method. 

To find the global error, we solve the recursion relation E k+ \ = ^kEk + L k , 
E Q = 0, cf. (27). The solution is 



(38) 



E k = 



-1 x^k-l 



E 



3=0 



Again, only the error in the nonstiff component propagates. 

Finally, we compute the error in the Evans function. Estimating the various 
terms in pip , we find that 

(39) E D =-§(« + «+)([£]i«+(0) + u(0) [E+} 2 ) + C(A _1 / 2 /i 8 , \- 3 ' 2 h 2 ). 

Hence, we wish to compute X = [E]iv+(0) + u(0) [E + ]2- Substitution of the exact 
solution, given in (fTfJ|) and (fT2)) . and the global error (j3"5)) yields 



X 



N 

X E 

3=0 



where 



N 



3=0 



3-1 



m + Lf + ) + K - 2 x 2 + o( K - 3 h A ). 



3-1 



X2 = H( lb 3- ha i E L * - L 3 $ +(°) + ^ + - ha t E L T + - L T + *- (o) 



i=0 



i=0 
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The sum Xj=o ( Z j + Z i + ) ^ s ^ e term within brackets in 132 j> . so again, it vanishes 
at all orders in h. For the <E>_ and <I> + terms, we use 

/< '• 1 rti+h i- 1 

(40) M60=/ M^«E/ ^-(x)dx = ^(L l a + /ia J ), 

i=0 i=0 

where the approximate equality becomes exact in the limit L — ► oo; again, we 
assume that L is so large that we can neglect any errors here. Substitution of (|4H|) , 
its analogue for $ + , and (|37a|) and (|37b|) yields 

N , N 

* 2 = -E ha i £ ° + + E K 1 " + 

j=0 ^ t=0 

+ ^ (hail? + hL?aj + ha+Lf + + hl a l ' + a+) J + 0(/i 8 ), 

where the remainder term comes from estimating terms like ^ • LfL?. This 
nested sum can be written as the product of two sums: 

N N 

X2 = E K«i + <4) ■ E ( z ? + z ^ + ) + °^ 8 )- 

i=0 3=0 

So we arrive again at the sum X^jLo ( Z j Z j + ) ' wm ch vanishes at all orders. 

Substituting everything back into (]39[) , we find that the error in the Evans func- 
tion is given by 

(41) E D = 0{\- l ' 2 h\\- l h\ \-^ 2 h 2 ). 

This is clearly better than the fourth-order Magnus method, with Ed = O(X h 4 ), 
and the exponential midpoint rule, with Ed = 0(X^ 1 ^ 2 h 2 ). 

5. Numerical experiments 

In this section, we evaluate the Evans function for a particular example. The 
error in this computation is determined and compared against the estimates derived 
in the previous sections. 

The example is the Fisher equation 

(42) U t = Mix +11- U 2 . 

This is a reaction-diffusion equation of the form (fT|). Fisher [T2] used it to de- 
scribe the transmission of genes in a population. It is now viewed as the prototype 
equation admitting travelling front solutions [35J §H-2]- 

The Fisher equation supports a travelling wave solution with wave speed c = 
— 1\/6- In fact, this solution is known analytically: 

u(x, t) = u(£) = j where £ = x — §\/f3t- 

(l + e^/V 5 ) 

Suppose that we wish to determine the stability of this travelling wave. We are led 
to consider the eigenvalue problem ([?]), which in this case reads 

1 ] 



(43a, | 
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where 
(43b) 



p(0 = i - 2&(0 = 1 - 



We solve this equation with the fourth-order Magnus method, given by (TT5)) . In the 



previous section, we derived the local error estimate 
this estimate evaluates to 



For the Fisher equation, 



(44) 



fe 5 e ;/V6(_ 8e 3 S /V6 + 33e 2;/y/6 + 702e ;/V6 + j) 

38880(1 + 6?/^) 6 
18(1 + e«/^) 3 



where K = y/ 'c 2 + 4(A + 1) and £ is short for 

The global error estimate is given in (|29|) . Assuming that £o is negative and so 
large in magnitude that we can take £o = — oo, we find that 

'V6h 4 e^^(36e 4 ^^ + 180e 3 «/^ + 364e 2 «/^ + 353e«/^ + l) 



(45) £ fe «~ 

K 



38880(1 + e«/v^) f 



18(1 + ei/Vey 

Finally, estimate p3|) for the error in the Evans function is 

V6 ,4 



(46) 



1080 



-0.002268 ft 4 



This estimate is independent of the parameter k. 

We perform some numerical experiments to check the validity of these estimates. 
First, we solve (|4"31 from £ = — 30 using the fourth-order Gauss-Legendre method 
with step size h = 0.02. We will refer to this solution as the "exact" solution. Then, 
we take the "exact" solution at £ = — 1 and do a single step with the fourth-order 
Magnus method with step size h = 0.2 or h = 0.1. The local error can now be 
determined by comparing the result of this single step against the "exact" solution; 
this local error is plotted in the top row of Figure [TJ together with the local error 
estimate (j44]). The horizontal axis in the plots shows the imaginary part of the 
eigenvalue parameter A, which varies from i to 10 s i in our experiments. 

The global error can be determined by solving (|4"3"]) from £ = —30 till £ = — 1 with 
the fourth-order Magnus method and comparing it against the "exact" solution. 
This results in the bottom row of Figure [TJ Finally, Figure [5] shows the difference 
between the Evans function as computed by the Magnus method and the "exact" 
value, compared against the estimate (|46|) . 

All graphs show that the error estimates agree well with the actual error when A 
is moderately large in magnitude. However, the numerical method starts to break 
down when |A| increases above 10 7 . 

The implementation used in the experiments is a straightforward Matlab code. 
One detail proved to be important, namely, the computation of the matrix ex- 
ponential in (|18p . The standard routine for this is called expm and uses Pade 
approximation combined with scaling and squaring. However, we found that an 
alternative approach based on the Schur decomposition and implemented in the 
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Figure 1. The solid lines in the graphs on the top row show the 
local error committed by the fourth-order Magnus method. The 
step size h is 0.1 and 0.2 for the line labelled 1 and 2, respectively. 
The dash lines show the local error estimate (pPT| . On the bottom 
row, the solid lines shows the global error and the dash lines show 
the estimate (|45|) . 

Matlab routine expmdemo3 works better in our case. Specifically, when using 
Pade approximation, the numerical method loses accuracy around |A| — 10 5 , as 
opposed to |A| = 10 7 for the Schur decomposition. Generally, the Schur decompo- 
sition runs into trouble when the matrix to be exponentiated is nearly defective, 
but in our case the eigenvalues are far apart, cf. (|22| . The reader is refered to the 
article by Moler and Van Loan [57J for an extensive discussion on this subject. 

When the experiment is repeated with the exponential midpoint rule (| 16[) and 
the fourth-order Gauss-Legendre method (|35p . the error in the Evans function is 
as plotted in Figure [31 We concluded in Section EH] that, because the exponential 
midpoint rule suffers less from order reduction, it is likely to have a smaller error 
than the fourth-order Magnus method if |A| is large. The error plots confirm this. 

For the fourth-order Gauss-Legendre method, which is the method that is more 
relevant in practice, the error in the Evans function is shown in the right half of Fig- 
ured We found the estimate ([4Tj) for the error, and the numerical results show that 
the term of order A _1 /i 4 dominates: the line 10 3 1 A | 1 /i. 4 tracks the graphs closely 
for |A| up to 10 4 (the coefficient 10~ 3 was not determined by any computation, 
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Error in the Evans function for h = 0.2 Error in the Evans function for h = 0.1 




Figure 2. The solid line shows the error in the Evans function 
as evaluated by the fourth-order Magnus method, while the dash 
line shows the error estimate f|46|) . The step size is h = 0.2 for the 
graph on the left and h = 0.1 for the graph on the right. 



Error in D(k) for exponential midpoint Error in D(X) for Gauss-Legendre 

10 4 , ■ 1 | : ■ 




Figure 3. The left graph shows the error in the Evans func- 
tion as evaluated by the exponential midpoint rule (|16p . while the 
right graph shows the same for the fourth-order Gauss-Legendre 
method |35|) . The step size is h = 0.2 for the curve labelled 2 and 
h = 0.1 for the curve labelled 1. The dotted lines in the second 
plot show 10- 3 h 4 /\X\. 



in contrast to the coefficient in (|46[) . but it was chosen to give a suitable match). 
When |A| > 10 4 , the error committed by the Gauss-Legendre method starts to 
increase erratically, following roughly the equation Ed = 10~ 12 y/\X\. This is likely 
due to round-off error, as \y(£,)\ is approximately v|A[. This suggests that the loss 
of accuracy in the fourth-order Magnus method is also due to round-off error, exac- 
erbated by ill-conditioning of the matrix exponential (compare with the influence 
of the method for computing the matrix exponential, as noted on[T5|. 
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6. Conclusions 

We found that the fourth-order Magnus method, when applied to the linear 
differential equation Q in the regime |A| ^S> l/h 2 , commits a global error of or- 
der h 2 /iy\\\ (relative to the exact solution). It is remarkable that the Magnus 
method converges at all. The convergence result [25] for the Magnus series men- 
tioned earlier guarantees convergence only when |A| < Tr/h 2 , so the usual conver- 
gence proof for the truncated series does not hold. However, the error analysis 
in this paper shows that the method does indeed converge for equations of the 
form (jll). 

Given that the method converges, it is remarkable that the order of the method 
drops. This is connected to the concept of stability. The Magnus method solves au- 
tonomous linear equations exactly. A fortiori, the numerical solution of Dahlquist's 
test equation y' = ay (with a € C) is stable if and only if the exact solution is sta- 
ble, meaning that it converges to as i — > oo. Hence, the Magnus method is 
A-stable and even L-stable (see, e.g., Hairer and Wanner [TH] for a definition of 
these terms). Nevertheless, the fourth-order Magnus method suffers from order 
reduction in the current setting, in which the equation is nearly autonomous. This 
may be connected to the fact that the fourth-order Magnus method is not B-stable. 
A simple counterexample is given by the equation y' — Ay with A(x) = [ "q 1 _\ ] 
for x < 3 and A(x) = for x > 3. This equation is contractive, but the 

numerical solution does not preserve contractivity as h increases above 6. In con- 
trast, the Gauss-Legendre method and the exponential midpoint rule are known to 
be B-stable, and they do not suffer from order reduction. We refer again to Hairer 
and Wanner [16] for a precise definition of B-stability and its connection to order 
reduction. 

Similar results were obtained by Hochbruck and Lubich |17j , who treated Magnus 
methods applied to semi-discretized Schrodinger equations. They could prove that 
the method converges even when there is no known convergence result for the 
untruncated Magnus series. Gonzalez, Ostcrmann and Thalhammcr [13 found that 
the exponential midpoint rule suffers from order reduction when applied to semi- 
discretized parabolic equations. The matrices in the semi-descretized equations 
considered by them have negative eigenvalues that are large in magnitude, just as 
the problem ((4]) treated here. 

However, the fourth-order Magnus method regains the full order when combining 
the solution of the differential equation |4} satisfying the boundary condition at 
£ = — oo with the one satisfying the condition at £ = +00 to form the Evans 
function. As is clearly shown both by the analysis and by the experiment, the error 
committed by the fourth-order Magnus is of order h A uniformly in A. Nevertheless, 
the Gauss-Legendre method is still superior: its error decreases as |A| increases. 

The same holds to a lesser degree for the exponential midpoint rule. The anal- 
ysis indicates that the error commited by this method is of order A -1 / 2 /) 2 . The 
numerical results for the exponential midpoint rule do not quite seem to agree with 
this, but they also show that the error decreases as a function of |A|; the reason for 
this discrepancy is unknown. Nevertheless, we can conclude that the second term 
in the Magnus expansion (|15p actually harms the numerical algorithm when A is 
large in magnitude. 

This suggests that the Right Correction Magnus Series, as proposed by Degani 
and Schiff [8], or the modified Magnus method, as proposed by Iserles [TI5], might 
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perform well on this problem. We ran some preliminary experiments with these 
methods, which showed that the error in the stiff component is greatly reduced and 
comparable to the error committed by the Gauss-Legendre method. However, the 
nonstiff component seems to suffer from round-off error. A full analysis of these 
methods warrants further investigation. 

As explained at the start of Section [21 our interest lies in the stability analysis 
for travelling waves for the reaction-diffusion equation ([1]). By energy estimates 
similar to those in Brin 5, §3.2], we find that the eigenvalues are contained in the 
wedge given by 

ReA< i C 2 + max 5 |/'(u(0)|, 
ReA+|lmA| < c 2 + max ? |/' («(£)) |. 

For the Fisher equation used as an example in the previous section, we have 
maxj = 1. Therefore, the analysis reported in this paper is of limited 

use when assessing the stability of the travelling wave for this equation. However, 
for other equations the wedge l|47p and the regime in which our analysis is valid 
may overlap. 

We mentioned in the introduction that this work was instigated by the paper 
of Aparicio, Malham and Oliver 3 . That paper studies a similar equation in the 
regime arg A = 7r, while the analysis in the current paper concerns the complemen- 
tary regime | arg A| < it — e with e > 0. In fact, it is possible to derive a combined 
estimate for the local error valid in both regimes using WKB-analysis. 

Finally, we wish to stress that the analysis reported here is only valid for scalar 
reaction-diffusion equations. We have not done a full analysis for systems of 
reaction-diffusion equations of the form ([1]). However, numerical experiments sug- 
gest that also in this case, the Magnus method suffers from order reduction but 
recovers the full order when matching the solutions to evaluate the Evans function. 
Generally, the Gauss-Legendre method still outperforms the Magnus method when 
A is away from the essential spectrum, but the difference is not so pronounced. 
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